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We calculate the condensate fraction and the condensate and non-condensate spatial and momen- 
tum distribution of the Bose-Hubbard model in a trap. From our results, it is evident that using 
approximate distributions can lead to erroneous experimental estimates of the condensate. Strong 
interactions cause the condensate to develop pedestal-like structures around the central peak that 
can be mistaken as non-condensate atoms. Near the transition temperature, the peak itself can 
include a significant non-condensate component. Using distributions generated from QMC simula- 
tions, experiments can map their measurements for higher accuracy in identifying phase transitions 
and temperature. 



The Bose-Hubbard (BH) model has been the focus of 
intensive research over the past decade as a prototypical 
example of strongly correlated physics, especially since 
the model was realized with cold atoms in optical lattices 
PP. Optical lattice experiments (OLE) are prime candi- 
dates for quantum simulations due to their extensive tun- 
ability and ease of control. Therefore they serve as ideal 
systems to study dynamical phenomena of many body 
effects in strongly interacting systems. Before studies 
can be meaningful, however, systematic characterization 
of equilibrium properties are crucial, the most important 
quantities being temperature (T) and density. However, 
direct or in-situ measurements of temperature in OLEs 
is an area of active research [2j |3]. Further, the study 
of phase transitions requires a good understanding of an 
observable that can be used as a probe for the state of 
the system. Theoretically, the natural choice is the or- 
der parameter, which for the BH model is the condensate 
fraction no = Nq/N where No is the total number of con- 
densed atoms and N the number of atoms. Alternatively, 
the superfluid fraction [4] [5] characterizes the transition, 
but this is not simple to measure in cold atom systems. 
The most easily accessible observables in experiments are 
the entropy and no that come from time-of-flight (TOF) 
measurements. The former is measured from TOF mea- 
surements while atoms are in the harmonic trap (without 
the lattice) and are then isentropically transferred into 
the lattice [3j |6j [7]. The latter could be measured di- 
rectly from TOF expansions after all fields are "snapped 
off." The no is particularly useful since, combined with 
entropy, it could be used for thermometry in experiments 

EE]. 

In homogeneous systems, no is given by a delta func- 
tion at the origin in momentum space. In TOF images, 
therefore, its signal would be visible from the appearance 
of a sharp peak. The presence of the parabolic trapping 
potential in OLEs, however, renders the system inhomo- 
geneous and no is no longer simply given by the occupa- 
tion number at zero momentum. 

The most common modelling approach to handle the 
trap is mean field (MF) theory: e.g. the Hartree-Fock- 



Bogoliubov-Popov approximation [9HTT] for small inter- 
action strengths (U) or the site-decoupled approach for 
large U [12], together with the local density approxima- 
tion (LDA). In situations where MF fails or quantitative 
comparison to measurement is important, we can resort 
to exact quantum Monte Carlo (QMC) techniques. QMC 
has been used to directly compare observables with mea- 
surements. However, to the best of our knowledge, the 
order parameter has not been computed and compared 
directly to experiment [T3] . 

In the most common approach, experimental TOF im- 
ages are heuristically fit to obtain the number of con- 
densed atoms under the peaks and the remaining non- 
condensed atoms. The ratio of the former to the sum 
of two is defined as the peak weight [8 or the coher- 
ence fraction [3 (/o) that serves as a proxy for no- In 
previous experiments, thermometry was done by compar- 
ing full momentum distributions together with /q, peak 
width (wq) and the visibility directly to QMC results 
[8] . The last three observables were further used to char- 
acterize the critical temperature (T c ) for the transition 
from the normal to the superfluid phase. Unfortunately, 
these probes are not necessarily reliable estimates of the 
order parameter since the relation between no and fo is 
not well understood. Previous comparisons show large 
differences [3]. This is unfortunate because no is a sim- 
ple probe for this transition and is also indicative of the 
effect of interactions, i.e. the quantum depletion (QD). 
Combined with the entropy measurements, no would be 
an excellent probe for the temperature. 

The difficulties in characterizing the mapping stems 
from the fact that the underlying condensate and non- 
condensate distributions are not well understood in 
trapped systems. If they were known, then no could 
be estimated by counting the number of atoms in the 
condensate peak and in the background. (We note that 
an analytical MFT method has been developed for a ho- 
mogeneous system [14]. Although qualitatively useful, 
this method finds limited use in real trapped systems 
where it must rely on LDA.) Another serious problem 
is the role of interaction during TOF expansion; some 
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have argued that the effect is small since by the time the 
wavefunctions from two adjacent lattice sites start over- 
lapping, the densities would have dropped dramatically 
[15] • However, others argue that the effects are signif- 
icant and lead to hydrodynamic effects in TOF images 
[3] . The crucial quantity that dictates the significance of 
these effects is the initial density. From two different sets 
of experimental data, it seems that for low densities (cen- 
tral filling of 1 or less) [8] interaction effects are small, 
while at large densities [3] (central filling of 3) they could 
be significant. Even in the case of non-interacting (NI) 
expansion, however, the interference pattern - reminis- 
cent of Fresnel diffraction - has to be explicitly accounted 
for [T5] . 

In this paper, we calculate no and the spatial and 
momentum distributions for several low density systems 
modelled with the BH Hamiltonian: 

h = -t + \ Yl ft< ( ft< ~ x ) ~ Yl w 
(ij) i i 

where t is the hopping integral between nearest neighbor 
sites i and j, a J (a^) is the Boson creation (annihilation) 
operator, hi = a\ai is the number operator, and U is the 
on-site repulsive interaction. Here, pi = p — V(r^/a) 2 
includes both the chemical potential term and spherical 
harmonic confining potential with a, the lattice spacing 
and V = \muj 2 the curvature that is given by the mass 
(m) and trap frequency (uu). Energies are given in atomic 
recoil energy units: E r ~ 167nK for 87 Rb and the laser 
wavelength A = 800 nm that is used to create the lattice. 
For our simulations, lattices are between 70 3 to 100 3 with 
open boundary conditions and the number of particles, 
N ~ 58000 to 64000. 

We calculate the single-particle density matrix 
Pi(hj) — (^1%) IH 111 HS1 HZ] using the stochastic se- 
ries expansion and the directed loop update algorithm 
[T8H20] . Then the occupation of the single particle states 
is defined by pi\ipi) — A^|^), where the largest eigen- 
value Nq of pi gives the number of condensed atoms and 
|^o) is the condensate wavefunction. The other atoms 
N nc — Ni = N — No are non-condensed atoms. 

For large systems, obtaining all the occupation modes 
is challenging because of the Monte Carlo noise in p\ and 
the complexity of a complete diagonalization of p\ . How- 
ever, since the condensate is not fragmented and occu- 
pies only one mode in these systems, we use an iterative 
diagonalization procedure to obtain the spatial conden- 
sate wave function, ^o( r ) an d the condensate momentum 
wave function, 0o(fc) = J-[ipo(r)] where T is the Fourier 
transform. For an extensive range of systems that we 
tested, the method is robust and is able to withstand 
statistical noise provided we do not enforce Hermitian 
symmetry; doing so biases the results. We use other sym- 
metries to reduce the noise in p\. For a full discussion, 
we refer to [21]. The spatial non-condensate distribution 



is then given as n nc (r) = n(r) — No\^o(r)\ 2 , where n(r) 
is the spatial density. The total momentum distribution 
is 

n(k) = \w(k)\ 2 Y,e iHj - l) p(i,j) =n (fc)+^> p (fc). 

jl p=l 

(2) 

We have explicitly included the wannier envelope \w(k)\ 2 
that is needed to go back to the continuum model 
from the lattice BH model. The last term on the rhs 
is the momentum non-condensate distribution n nc (k). 
In order to match with experiments, we must include 
the finite TOF effects [15] , by adding an additional 
site dependent phase term to ([2| so that nj(fc) = 

|w(fc)| 2 Wp£ j7 e <fc ^^ 

r is the TOF time, m is the particle mass. We use rio(fc) 
to denote the finite TOF condensate and n^ c (k) for the 
finite TOF non-condensate distributions. 

In fig. ([lji-d), we present the different density pro- 
files corresponding to large U. (The trap frequency has 
been adjusted so that n(r = 0) = 1; the density is 1 
atom/lattice site at the center of the trap.) We note that 
the underlying distributions here are exact. At U/t = 40, 
the Mott Insulating (MI) domain appears as an integer 
plateau. Similar distributions obtained in previous stud- 
ies [9HT2] using MFT + LDA cannot capture the inter- 
face between the developing MI and superfluid domains 
that play an increasing role at high U/t since these stud- 
ies effectively stitch together densities for different radii, 
while assuming that there are no exchanges across differ- 
ent shells. Therefore, studies using these distributions as 
masks for fitting purposes must be carefully interpreted. 

The column integrated momentum distributions in the 
far-field, shown in fig. ([]J-h), exhibit secondary peaks 
similar to [22j [23]. Note that these peaks can arise 
not only due to the finite extent of the condensate but 
rather all modes (fig. ([]J-h) insets). However, it is only 
around the MI regime that it will be primarily due to the 
condensate, and so we see that the peak forms around 
k — 27r/£o, where £0 is the width of the condensate. This 
is specifically due to the way the condensate forms be- 
tween boundaries (in this case between the MI domain 
and the vacuum). We present statistics of the secondary 
peak in table Q from which it is evident that its size 
relative to the central peak and the contribution of the 
condensate mode to it is much larger (2.8 times at low T 
and 5 times at high T) when the system has a MI domain. 

Finite TOF effects, presented in fig. ([l}-l), alter the far 
field distributions by suppressing and blurring the cen- 
tral low k values, as expected [15] . Higher order modes 
of pi with rapid spatial variations are not significantly 
affected by the site dependent phase shift. Thus, the 
maximum effect is on the condensate distribution. The 
time scale for the condensate to reach the far field (jff) 
is oc i?£o(l — £o/2-R), where R is the radial extent of 
the condensate. This leads to larger Tff for U/t = 25 
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U/t = 25 N ~ 60, 000 uj = 68.1Hz 



U/t = 40 N ~ 58, 000 u = 67.6Hz 



k b T/t = 2.456 n = 0.218(2) 
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FIG. 1. Top row (a-d) are the spatial particle densities, middle row (e-h) and bottom row (i-1) are the far field column 
integrated and finite TOF momentum distributions. Insets on the middle line (e-h) show fine features for all distributions and 
last line (i-1) show zoomed versions of the finite TOF non-condensate distribution (n^ c (fc)). For the momentum distributions 
in the middle and last line the y-axis is in arbitrary units, but the scale is the same for all images. 



and so for the fixed r = 20ms, the central peak sees a 
greater suppression (and surrounding region greater en- 
hancement) than U/t = 40. Although, n(k) and n T (k) 
are both broader due to the secondary peaks, the latter 
case has relatively more condensate atoms. They would 
not be captured in fitting schemes used in experiments 

EE]. 

We note that the broader structure is observable 
within the experimental resolution. Using Afc = 
(m\7r/hr)Ar/a^ where A the wavelength of the opti- 
cal lattice, r the expansion time, Ar the resolution, 
we obtain the Afc resolution. Using Rb 87 , r = 20 ms, 
A = 800nra and typical resolving power of Ar = 3/ira 
gives Afca ~ 0.0267r/a. Features in fig. Q-l) are spread 
over Afc ~ 0.067r/a. 

Fig. ([2^ and b) show underlying distributions around 
the transition temperature (T c ) for a weakly interact- 
ing (U/t = 3.4) and a strongly interacting (U/t = 25) 
system when n nc = 0.98(5) and 0.998(6). In accordance 
with previous studies, we see strong central peaks [22) 124] 
even at these warm temperatures. Furthermore, the cen- 
tral and surrounding Bragg-peaks are not representative 
of the condensate on its own. The magnitude of the error 
due to approximate fits is specific to the fitting procedure 
used. However, we can systematically analyze boundaries 
of possible errors due to the issues we discussed above. 



We estimate the coherence fraction by calculating the ra- 
tio of atoms under the central peak (up to a limit ki) to 
the total number of atoms. (The numbers are calculated 
from column integrated images without the wannier en- 
velope so that we do not worry about its subtraction.) 
The range is chosen large enough to accommodate finite 
optical resolvability in experiments. It should also allow 
for related fitting procedures. 



U/t = 3.4 u = 37.7 Hz 



U/t = 25 uj = 68.1 Hz 



^ xlOl 



k b T/t = 3.93 n = 0.0151(1) k b T/t = 3.27 n = 0.00142(4) 





FIG. 2. Closer look at the distributions near the transition 
temperature for a system of N ~ 64, 000 atoms, (a) and (b) 
are finite time TOF distributions with r = 20ms. Depending 
on the cut-off k value chosen (ki) there are between 500 to 
5000 atoms for U/t = 3.4 (ki = 0.01,0.1) and 300 to 1260 
atoms for U/t = 25 (ki = 0.035,0.075). The exact number is 
960 and 90 respectively. The insets show n nc (k) and no(fc) 
with 3 peaks. 
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FIG. 3. Condensate fraction (no) as a function of temperature 
(T) compared with coherence fraction (fo(T)) measurements 
for different cutoff fa. (fo = 5^|fc|<fc no(k)/N.) We have 
fit the no to the functions given by ^ where Uq = 0.638(7) 
(0.502(6)), k b T t /t = 2.99(2) (2.36(1)) and g = 1.91(9) (2.4(1)) 
for U/t = 25 (40). The short arrows indicate fa,T c /£ = 3.023 
(2.285) for U/t = 25 (40) where n - 0.01 

We present comparisons with exact results in fig. ([3^l 
and b) approximated by: 

n (T) = nS(l - exp(g(l - ^)), (3) 

away from the critical regime [25]. We see that at low T, 
the depletion is severely overestimated whereas around 
the transition it is underestimated. It is worth noting 
that general trend of the error is in accordance with fo 
measurements presented in [3] , where it is further exacer- 
bated by experimental noise and additional fitting errors 
due to inaccurate background subtraction and possible 
interaction effects during expansion of the gas. 

Care is needed in the estimation of T c using fo. In 
the general case, extrapolating from our data, for ki > 
0.035, fo(h) > no. If T c is to be estimated correctly, 
\d 2 fo/dT 2 \ > \d 2 no/dT 2 \ is required for rapid conver- 
gence to the same uq(T) — >> 0. Futher, following [26], a 
study of the trap-size scaling behavior analogous to fi- 
nite size scaling studies done for homoegeneous systems 
would be needed to exactly identify critical exponents 
and T c for trapped systems. Here we use a simpler 
working definition no(T c ) ~ 0.01 that suffices to show 
that estimates of T c using fo would be erroneous: in 
fig. pk and b), k b T c /t = 3.022 (2.284), but using / , 



3.023(2.29) < k b T c /t < 3.3(2.81) for U/t = 25 (40). Fur- 
thermore, it should also be clear that a smooth change 
in fo(T) across the transition might make it impossible 
to use the inflection point as an accurate estimate of T c . 

A simple way for experiments to obtain the correct fit- 
ting function might be to calibrate the coldest possible 
measurements of fo against Uq ~ no(T = 0), the satu- 
rated no from QMC. The approximate temperature could 
be approximated from entropy measurements as in |3j [7] . 
The TOF image used to estimate fo should be then fitted 
correctly: the bi-modal structure of the condensate seems 
amenable to be fit by a narrow Thomas-Fermi like func- 
tion plus a broader Gaussian. Fitting the non-condensate 
correctly is somewhat easier since, at low T the height 
of n nc (k) at low k is small and fluctuations can be inte- 
grated out. On the other hand, at high T it is flat under 
the peak and cannot be assumed to be gaussian-like out- 
side. 

We have presented the effects of interactions on the 
components of the spatial and momentum distributions 
of bosonic particles in a trapped optical lattice. Our un- 
biased estimates of pi elucidate the potential problems 
in mapping between the coherence fraction and the exact 
no, and how to account for them. Using exact distribu- 
tions for QMC, experiments will get more accurate esti- 
mates of no- In a future work, we will study the effects 
of interaction during TOF and present comparisons with 
experimental systems at higher density [27] . We will also 
explore the entropy (S) mapping using QMC that is need 
to compare with no(S) from experiments. The method 
we have presented in this paper will be crucial to such 
studies. 
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TABLES 



U/t 


k b T/t 


n (kp2)/n(kp 2 ) 


no(k p2 )/n(0) 


25 


2.46 


0.405 


0.0098 


0.98 


0.742 


0.0127 


40 


1.96 


0.696 


0.0414 


0.98 


0.86 


0.034 



TABLE I. Statistics of the secondary peaks: fraction of con- 
densate at the maxima of the secondary peak (k P 2) and size 
relative to the central peak. 



